Spatio-temporal correlations can drastically change the response of a MAPK pathway 
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Multisite covalent modification of proteins is omnipresent in eukaryotic cells. A well-known ex- 
ample is the mitogen-activated protein kinase (MAPK) cascade, where in each layer of the cascade 
a protein is phosphorylated at two sites. It has long been known that the response of a MAPK 
pathway strongly depends on whether the enzymes that modify the protein act processively or dis- 
tributively: a distributive mechanism, in which the enzyme molecules have to release the substrate 
molecules in between the modification of the two sites, can generate an ultrasensitive response and 
lead to hysteresis and bistability. We study by Green's Function Reaction Dynamics, a stochastic 
scheme that makes it possible to simulate biochemical networks at the particle level and in time and 
space, a dual phosphorylation cycle in which the enzymes act according to a distributive mechanism. 
We find that the response of this network can differ dramatically from that predicted by a mean- 
field analysis based on the chemical rate equations. In particular, rapid rebindings of the enzyme 
molecules to the substrate molecules after modification of the first site can markedly speed up the 
response, and lead to loss of ultrasensitivity and bistability. In essence, rapid enzyme-substrate 
rebindings can turn a distributive mechanism into a processive mechanism. We argue that slow 
ADP release by the enzymes can protect the system against these rapid rebindings, thus enabling 
ultrasensitivity and bistability. 
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I. INTRODUCTION 



Mitogen-activated-protein kinase (MAPK) cascades 
are ubiquitous in eukaryotic cells. They are involved 
in cell differentiation, cell proliferation, and apoptosis 
Q. MAPK pathways exhibit very rich dynamics. It 
has been predicted mathematically and shown exper- 
imentally that they can generate an ultrasensitive re- 
sponse 0, H, 01 and exhibit bistability via positive feed- 
back 0]. It has also been predicted that they can gen- 
erate oscillations [1, 0) Hi, amplify weak but attenuate 
strong signals [§], and give rise to bistability due to en- 
zyme sequestration [id MAPK pathways are in- 
deed important for cell signalling, and for this reason 
they have been studied extensively, both theoretically 

i i i a i i [MM 13, JMIJIJ^ m, m 

and experimentally 0, 0, [a, 0, [M, [20, [211, [23 • However, 
in most theoretical analyses, the pathway is modelled us- 
ing chemical rate equations 0, i, [1, [1 [Tfl [H, [H [H, [ll] . 

This is a mean-field description, in which it is assumed 
that the system is well-stirred and that fluctuations can 
be neglected. Here, we perform particle-based simula- 
tions of one layer of the MAPK cascade using our re- 
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cently developed Green's Function Reaction Dynamics 
algorithm [IE [IS- Our simulations reveal that spatio- 
temporal correlations between the enzyme and substrate 
molecules, which are ignored in the commonly employed 
mean-field analyses, can have a dramatic effect on the 
nature of the response. They can not only speed up 
the response, but also lead to loss of ultrasensitivity and 
bistability. 

The response time, the sharpness of the input-output 
relation, and bistability are key functional characteris- 
tics of signal transduction pathways. The response time 
does not only determine how fast a cell can respond to 
a changing environment, but has also been implicated 
to underlie many cellular decisions. For example, pro- 
cesses such as cell proliferation and differentiation, se- 
lection of T cells, apoptosis, and cell cycle progression 
are believed to be regulated by the duration of the sig- 
nal [H, [H [IS [H, [13, [H, [li- The sharpness of the 
input-output relation, or the gain, is a key property of 
any signal transduction pathway, since it directly affects 
the signal-to-noise ratio. Bistability can lead to a very 
sharp, all-or-none response [5,], buffer the cell against fiuc- 
tuations in an input signal, and makes it possible to lock 
the cell in a given state. Indeed, bistability, or more in 
general multistability, plays a central role in cell differ- 
entiation [s^, [3l| . It is thus important to understand the 
mechanisms that underlie bistability, the gain and the 
response time of MAPK pathways. 

A MAPK cascade consists of three layers, where in 
each layer a kinase activates the kinase of the next layer. 
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Importantly, full activation of the kinase requires that 
it becomes doubly phosphorylated (see Fig. [T|). Kinase 
activation is regulated via a dual phosphorylation cycle, 
in which the upstream kinase and a phosphatase control 
the phosphorylation state of the two sites of the kinase in 
an antagonistic manner. A key question is whether the 
enzymes that modify the kinase act in a processive or in 
a distributive manner [2,, i4] ■ In a distributive mecha- 
nism, the enzyme has to release the substrate after it has 
modified the first site, before it can rebind and modify the 
second site. In contrast, in a processive mechanism, the 
enzyme remains bound to the substrate in between the 
modification of the two sites. While a processive mech- 
anism requires only a single enzyme-substrate encounter 
for the modification of both sites, a distributive mecha- 
nism requires at least two enzyme-substrate encounters. 

Mean-field analyses based on the chemical rate equa- 
tions have revealed that whether the enzymes act ac- 
cording to a processive or a distributive mechanism has 
important functional consequences for the response of 
a MAPK pathway. A distributive mechanism can gen- 
erate an ultrasensitive response since the concentration 
of the fully activated kinase depends quadratically on 
the upstream kinase concentration 0, y, Moreover, 
if the enzymes are present in limiting amounts, enzyme 
sequestration can lead to bistable behavior if they act 
distributively These mean- field analyses, however, 

assume that at each instant the molecules are uniformly 
distributed in space. Here, we show using particle-based 
simulations that spatio-temporal correlations between 
the enzyme and the substrate molecules can strongly af- 
fect the response of a MAPK pathway. 

We perform particle-based simulations of one layer of 
a MAPK pathway in which the enzymes act according to 
a distributive mechanism. The simulations reveal that 
after an enzyme molecule has dissociated from a sub- 
strate molecule upon phosphorylation of the first site, 
it can rebind to the same substrate molecule to modify 
its second site before another enzyme molecule binds to 
it. Importantly, the probability per unit amount of time 
that such a rebinding event occurs does not depend upon 
the enzyme concentration. As a result, enzyme-substrate 
rebindings can effectively turn a distributive mechanism 
into a processive one, even though modification of both 
sites of a substrate molecule involves at least two colli- 
sions with an enzyme molecule. Indeed, a distributive 
mechanism not only requires a two-collision mechanism, 
it also requires that the rates at which they occur depend 
upon the concentration. 

These rebindings have important functional conse- 
quences. Since rebindings effectively turn a distribu- 
tive mechanism into a processive one, ultrasensitivity and 
bistability via enzyme sequestration are lost. Moreover, 
rebindings strongly reduce the gain of the network. We 
investigate in depth the scenarios in which rebindings 
become important. This reveals that the importance of 
rebindings depends on the concentration and the diffu- 
sion constant of the molecules: the lower the concentra- 



tion and/or the diffusion constant, the more likely an 
enzyme molecule rebinds a substrate molecule to mod- 
ify the second site before another enzyme molecule does. 
Since enzyme-substrate rebindings are faster than ran- 
dom enzyme-substrate encounters, this observation leads 
to the counter-intuitive prediction that slower diffusion 
can lead to a faster response. We also find that the im- 
pact of rebindings strongly depends on the time it takes 
to re-activate the enzyme after it has modified the first 
site. If, for instance, the ADP/ATP exchange on a ki- 
nase has to take place after the kinase has dissociated 
from the substrate upon phosphorylation of the first site, 
but before it can bind the substrate again to modify the 
second site, then either slow ADP release or slow ATP 
supply will make enzyme-substrate rebindings less im- 
portant. ADP release from protein kinases has been re- 
ported to be fairly slow [32], suggesting that slow ADP 
release might be critical for generating ultrasensitivity 
and bistability. 

The importance of rebindings relies on the interplay 
between reaction and diffusion at short and long length 
and time scales. This means that the algorithm should 
correctly capture the spatio-temporal dynamics of the 
system at both scales. In this manuscript, we present 
and apply an enhanced version of our recently developed 
Green's Function Reaction Dynamics algorithm. This 
particle-based algorithm is not only even more efficient 
than the original GFRD scheme, which is already 4 to 5 
orders more efficient than brute-force Brownian Dynam- 
ics 23], it is also exact. 

Biological systems that exhibit macroscopic concentra- 
tion gradients or spatio-temporal oscillations, which have 
recently been studied extensively, are typically consid- 
ered to be reaction-diffusion problems. We believe that 
our simulations are the first to show that in a biological 
system that is spatially uniform, spatio-temporal corre- 
lations on molecular length scales can drastically change 
the macroscopic behaviour of the system. This under- 
scores the importance of particle-based modelling of bio- 
logical systems in time and space. 



II. MODEL 

A. Dual phosphorylation cycle 

We consider one layer of the MAPK pathway, con- 
sisting of one dual modification cycle, as shown in Fig. 
[T] Phosphorylation and dephosphorylation proceed via 
Michaelis-Menten kinetics and according to an ordered, 
distributive mechanism. Importantly, we assume that 
the enzymes are inactive after they have released their 
modified substrate; before they can catalyse the next re- 
action, they first have to relax back to the active state. 
The inactive state could reflect that the enzyme is in 
an inactive conformational state after it has released its 
product. For the kinase it could also reflect that after 
it has released its substrate, ADP is bound; only when 
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FIG. 1: Dual phosphorylation cycle of one layer of the MAPK 
cascade. MAPK (K) is activated via double phosphorylation 
by the kinase MAPKK (KK) of the upstream layer and de- 
activated via dephosphorylation by a phosphatase (P). It is 
assumed that the enzymes KK and P act distributively and 
become inactive (KK* and P*) immediately after the sub- 
strate has been modified, releixing back to the active state 
with a characteristic time scale Trei. 



ADP has been released and ATP has been bound, does 
the enzyme become active again. As we will discuss in 
detail below, the timescale for re-activation, Tioi, plays a 
key role in the dynamics of the system. 

This model is described by the following reactions: 
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The first two reactions describe the phosphorylation of 
the kinase of interest, MAPK (K), by the upstream ki- 
nase, MAPKK (KK), while Eqs. [3] and H describe its 
dephosphorylation by the phosphatase (P). The inactive 
state of the enzymes after they have released their prod- 
uct is denoted by the superscript *, and the relaxation 
towards the active state is described by the last two equa- 
tions. For simplicity, we assume that re-activation can be 
described as a simple unimolecular reaction with a time 
scale Trci — l/fc?. We also assume that the system is 
symmetric, meaning that the rate constants for the phos- 
phorylation reactions are equal to the corresponding rate 
constants for the dephosphorylation reactions. We will 
systematically vary the relaxation time Troi, and the con- 
centration and the diffusion constant of the particles, D 
(see below). For the other parameter values, we have 
taken typical values from the literature (see Methods). 



B. Green's Function Reaction Dynamics 

We will compare the predictions of a mean-field model 
based on the chemical rate equations [l0[ with those of 
a model in which the particles are explicitly described in 
time and space. In this particle-based model, it is as- 
sumed that the molecules are spherical in shape, have 
a diameter a, and move by diffusion with a diffusion 
constant D. Moreover, two reaction partners can react 
with each other with an intrinsic rate = ki or ki^, re- 
spectively, once they are in contact, and two associated 
species can dissociate with an intrinsic dissociation rate 
kd = fc2 or ^5, respectively. 

One algorithm to simulate this particle-based model 
would be Brownian Dynamics. However, since the con- 
centrations are fairly low, much CPU time would be 
wasted on propagating the reactants towards one an- 
other. We therefore employ our recently developed 
Green's Function Reaction Dynamics algorithm, which 
uses Green's functions to concatenate the propagation of 
the particles in space with the chemical reactions between 
them, allowing for an event-driven algorithm [23, . ,24,] (see 
Methods). 

III. RESULTS 

A. Rebindings 

To understand the response of the dual phosphoryla- 
tion cycle, it is critical to consider the distribution of 
association times for a bimolecular reaction. We con- 
sider a simple bimolecular reaction, A -\- B ^ C , where 
one A molecule can react with one of TV B molecules to 
form a C molecule in a volume V . A model in which it 
is assumed that the particles are uniformly distributed 
in space at all times, be it a mean-field continuum or a 
stochastic discrete model, predicts that this distribution 
is exponential (see Fig. [2]). In contrast, in a spatially- 
resolved model, the distribution of association times is 
algebraic on short times and exponential only at later 
times 33]. 

The difference between the well-stirred model and the 
spatially-resolved model is due to rebindings. In a well- 
stirred model, the propensity that after a dissociation 
event the A molecule reacts with a B molecule only de- 
pends on the total density N/V of B molecules, and 
not on their positions — in a spatially resolved model this 
would amount to putting the dissociated B particle to a 
random position in the cell. Since the total density of 
B is constant, the association propensity is constant in 
time, leading to an exponential waiting-time distribution 
in the well-stirred model. In the spatially resolved model 
the situation is markedly different. The B molecule that 
has just dissociated from the A molecule is in close prox- 
imity to the A molecule. As a consequence, it can rapidly 
rebind to the A molecule before it diffuses away from it 
into the bulk. Such rebindings lead to the algebraic de- 
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FIG. 2: The distribution of association times for a bimolec- 
ular reaction for different concentrations. The system con- 
sists of one A molecule that can associate with A'' = 10 
B molecules according to the reaction A + B ^ C. For 

t < Tmoi ~ cr^ / D ~ Ifis the distribution decays as t~^^'^, 
for Tnioi < t < Tbuik ~ Is it decays as t~^^^ , while for t > rbuik 
the distribution decays exponentially. The algebraic decay is 
due to rebinding events, in which a dissociated B molecule 
rebinds the A molecule before diffusing into the bulk; this 
is unaffected by the concentration of B, [B]. The exponen- 
tial relaxation is due to B molecules that arrive at A from 
the bulk, and is a function of [B] . The concentration was con- 
trolled by changing the volume from V = 1, 0.33, and 0.1 /^m^, 
corresponding to [B] = N/V = 16, 48, and 160 uM, respec- 
tively, fca = 0.056 nM-^s-^ fed = 1.73 s'^ D = l/im^s^^ and 
(T = 5 nm. 



FIG. 3: Average response time as a function of the diffusion 
constant D for t^^i ~ 10 ms (blue line) Troi = 1 ^s (green line), 
as predicted by the particle-based model; for comparison, the 
predictions of the mean-field model based on the ODE chem- 
ical rate equations are also shown (black lines). Initially, only 
the phosphatases are active; at t — the upstream kinases 
are activated, and plotted is the time it takes on average to 
reach 50% of the final steady-state level of doubly phospho- 
rylated substrate (Kpp). The optimum in the particle-based 
model is due to the interplay between phosphorylation of the 
first site, which slows down with decreasing diffusion constant 
since enzyme and substrate have to find each other at random, 
and phosphorylation of the second site, which speeds up with 
decreasing diffusion constant, because of enzyme-substrate re- 
bindings. 



B. Rebindings can speed up the response 



cay of the association-time distribution at short times. 
For times shorter than the time to travel a molecular 
diameter, t < t,„o1 ~ cr^/D (see Supporting Informa- 
tion), the dissociated B particle essentially experiences a 
surface of the A particle that is fiat, and its rebinding 
dynamics is given by that of a ID random walker return- 
ing to the origin, leading to the t~^^^ decay. At times 
Tmoi < t < Tbuik, the dissociated B particle sees the en- 
tire sphere of A, and the probability of a re-encounter 
event is that of a 3D random walker returning to the ori- 
gin, decaying as t~^^'^. At times t > Tbuik, the dissociated 
B particle has diffused into the bulk, and it has lost all 
memory where it came from. The probability that this 
molecule, or more likely, another B molecule binds the A 
molecule, now becomes constant in time, leading to an 
exponential waiting-time distribution at long times [s^ . 

Fig. [2] shows that the association-time distribution de- 
pends on the concentration for t > Tbuik, but not for 
t < Tbuik- Indeed, while the encounter rate between two 
molecules in the bulk depends on their concentration, the 
rate at which a rebinding event occurs is independent of 
it. As we will show below, this has major functional con- 
sequences for the response of the dual phosphorylation 
cycle. 



Fig. [3]shows the average response time as a function of 
the diffusion constant, for two different values of the life- 
time of the inactive state of the enzymes, Trei. The figure 
reveals that both the mean-field (ODE) and the particle- 
based model predict that there is an optimal diffusion 
constant that minimizes the response time. However, in 
the mean-field model the optimum is barely noticeable 
{47] [34] . To a good approximation, the mean- field model 
predicts that the response time increases with decreasing 
diffusion constant, because enzyme-substrate association 
slows down as diffusion becomes slower. In contrast, the 
particle-based model shows a marked optimum, which is 
most pronounced when Tj-ei is short. Clearly, the particle- 
based simulations predict that slower diffusion can lead 
to a faster response. 

The speed up of the response with slower diffusion is 
due to the interplay between enzyme-substrate rebind- 
ings, and enzyme re-activation. This interplay mani- 
fests itself in the distribution of the second-association 
time, defined as the time it takes for a substrate molecule 
that has just been phosphorylated (Kp) to bind a kinase 
molecule (KK) for the phosphorylation of the second site 
(Fig. [1]). After a kinase molecule (KK) has phospho- 
rylated the first site of a substrate molecule (K), it will 
dissociate from it. After dissociation, it is still in close 
proximity to the substrate molecule, and it will there- 
fore rapidly re-encounter the substrate molecule before 
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it diffuses away into the bulk. When the hfetime Trci of 
the inactive state of the kinase molecule is short com- 
pared to the time t^o\ it takes for the enzyme and sub- 
strate molecule to diffuse away from each other, the prob- 
ability that upon a re-encounter the enzyme molecule 
has become active again such that it can actually re- 
bind the substrate molecule, will be large. Hence, when 
Trci < Tmoi, the kinasc will often rapidly rebind the sub- 
strate molecule, leading to the characteristic algebraic 
decay of t"^/^ for Tmoi < t < rbuik (Fig. SK). How- 
ever, there is also a probability that the enzyme molecule 
will escape into the bulk before it rebinds the substrate 
molecule. If this happens, most likely another kinase 
molecule binds the substrate molecule. This scenario un- 
derlies the exponential form of the second-association- 
time distribution at longer times, with the corner time 
Tbuik ~ 0.1 — 10s. It can now also be understood why the 
marked peak in the distribution at short times (Fig. |4lA_) 
disappears when the enzymes' reactivation time Trci be- 
comes significantly longer than Tmoi (Fig. HP): after phos- 
phorylation of the first site, the kinase will rapidly re- 
encounter the substrate molecule many times, but since 
the enzyme is most probably still inactive, it cannot re- 
bind the substrate molecule, and it will therefore dif- 
fuse into the bulk. In the Supporting Information we 
derive analytical expressions for the enzyme-substrate 
rebinding-time distributions, and elucidate the different 
scaling regimes that can be observed. 

To understand why slower diffusion can lead to a faster 
response when the lifetime of the enzymes' inactive state 
is short (Fig. [3]), it is instructive to consider how the 
distribution of second-association times depends on the 
diffusion constant. Fig. [4lA_ shows that the corner at 
Tbuik shifts to longer times as the diffusion constant is 
decreased. This is because the rate at which a kinase 
molecule from the bulk encounters a given substrate 
molecule is given by l/rbuik = fco = 47ro-(I?E + £'s)[KK], 
where a is the sum of the radii of the enzyme and 
substrate molecules and De and Ds are the diffusion 
constants of the enzyme and substrate molecules, re- 
spectively. Clearly, substrate phosphorylation by kinase 
molecules that have to find the substrate molecules at 
random slows down as the molecules move slower. How- 
ever, the figure also shows that the distribution at the 
corner time of Tbuik decreases in magnitude while the 
peak at Tmoi increases in magnitude when diffusion be- 
comes slower. This means that as the diffusion con- 
stant becomes lower, phosphorylation of the second site 
is increasingly dominated by enzyme-substrate rebind- 
ings rather than by random enzyme-substrate encoun- 
ters. The probability that the enzyme molecule is still 
in the vicinity of the substrate molecule after it has re- 
laxed back to the active state, increases as the diffu- 
sion constant decreases, making a substrate-rebinding 
event more likely. This is demonstrated quantitatively 
in Fig. HP, which shows the probability that both sites 
on the substrate are phosphorylated by the same ki- 
nase molecule. As expected, this probability not only 



increases with decreasing lifetime of the enzymes' inac- 
tive state, but also with decreasing diffusion constant. 
Since enzyme-substrate rebindings are more rapid than 
random enzyme-substrate encounters, this explains why 
slower diffusion can lead to a faster response. 

While slower diffusion speeds up the modification of 
the second site by making rapid enzyme-substrate re- 
bindings more likely, it also slows down the modification 
rate of the first site since that is determined by the rate 
at which enzyme molecules find the substrate molecules 
from the bulk. This is the origin of the optimum diffusion 
constant that minimizes the response time (Fig. [3]). 



C. Enzyme-substrate rebindings can weaken the 
sharpness of the response 

Fig. E] shows the effect of enzyme-substrate rebind- 
ings on the steady-state input-output relation. It is seen 
that when the re-activation time of the enzymes is long, 
Tioi = 10 ms, the input-output relation is strongly sig- 
moidal (Fig. 03). Moreover, it does not depend much 
on the diffusion constant of the molecules, and it agrees 
quite well with that predicted by the mean-field model 
based on the chemical rate equations (Fig. [5j3 In con- 
trast, when Troi is short, i.e. Trei — 1/is, the input-output 
relation markedly depends on the diffusion constant (Fig. 
[5j^). For large diffusion constants, the response curve 
agrees well with that predicted by the mean-field model 
of a distributive mechanism. But for lower diffusion con- 
stants, it increasingly deviates from the mean-field pre- 
diction, and it becomes significantly less sigmoidal. 

It is commonly believed that multi-site covalent mod- 
ification can lead to a sigmoidal, cooperative response 
when the enzymes act distributively, but not when they 
act processively 0, [s^ . While in a distributive scheme 
modification of n sites of a substrate molecule requires 
at least n enzyme-substrate binding events, in a proces- 
sive scheme only one enzyme-substrate binding event is 
needed. This is often presented as the explanation for 
why a distributive mechanism enhances the sensitivity of 
the modification level to changes in enzyme concentra- 
tion. However, Fig. [SJ^ shows that when the enzymes' 
re-activation time is short and the species' diffusion con- 
stant is low, the input-output relation of a distributive, 
dual phosphorylation cycle approaches that of a proces- 
sive, dual phosphorylation cycle. This is due to enzyme- 
substrate rebindings. Even though during a rebinding 
trajectory the enzyme molecule is detached from the sub- 
strate molecule and two binding events are required for 
full substrate modification, the rate at which the second 
site is modified does not depend on the enzyme concen- 
tration (Fig. m . The sharpness of the response increases 
with the number of required enzyme-substrate binding 
events, but only when these depend on the enzyme con- 
centration. Enzyme-substrate rebindings effectively turn 
a distributive mechanism into a processive mechanism. 
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FIG. 4; Distribution of times it takes for a substrate that has just been phosphorylated once (Kp) to bind a kinase molecule 
(KK), for different diffusion constants. The enzyme reactivation time is (A) Trd ~ 1/is or (B) Trd ~ fOms; in both cases 
Tmoi ~ l^s and Tbuik ~ O.f — lOs. (C) Probability that the two modification sites of a substrate molecule are phosphorylated 
by the same kinase molecule as a function of diffusion constant, for different enzyme re-activation times Trei. It is seen that the 
probability of an enzyme-substrate rebinding event increases not only with decreasing enzyme reactivation time, but also with 
decreasing diffusion constant. The latter explains why slower diffusion can lead to a faster response, as seen in Fig. [3] 
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FIG. 5: Steady-state input-output relations for different dif- 
fusion constants when (A) Trei = 1 ^ls and (B) Trei = 10 ms. 
For comparison, we also show the predictions of a mean-field 
(ODE) model of a distributive system with D — 1 /im^/s, and 
that of a particle-based model of a processive system with 
D = 0.06/im^/s (only in panel A). Note that when the re- 
activation time Trci is short, the input-output relation of a 
distributive system approaches that of a processive one as 
the difi^usion constant is lowered (panel A; blue lines). 



vation time t^cI is long, spatio-temporal correlations are 
not important, and the system indeed exhibits bistability. 
But when is short, the probability that a substrate 
molecule that has just been phosphorylated once will be 
phosphorylated twice is larger than that it will be de- 
phosphorylated again: the chance that it will rebind the 
kinase molecule that has just phosphorylated it, will, be- 
cause of the close proximity of that kinase molecule, be 
larger than the probability that it will bind a phosphatase 
molecule, even though in this state there are many more 
phosphatase than kinase molecules to which the substrate 
molecule could bind to. These rebindings, or more pre- 
cisely, spatio-temporal correlations between the enzyme 
and substrate molecules, are the origin of the loss of bista- 
bility when Tici is short (Fig. 



The effect of concentration 



D. Rebindings can lead to loss of bistability 

Markevich et al. have shown that bistability can 
arise in a dual phosphorylation cycle when the enzymes 
act distributively and are present in limiting concentra- 
tions [l^. The idea is that if the substrate molecules 
are, for example, predominantly unphosphorylated and 
a substrate molecule is phosphorylated to become singly 
phosphorylated, it will most likely bind a phosphatase 
molecule to become unphosphorylated again, instead 
of a kinase molecule to become fully phosphorylated — 
when most of the substrate molecules are unphosphory- 
lated, the kinase molecules are mostly sequestered by the 
unphosphorylated substrate molecules, while the phos- 
phatase molecules are predominantly unbound. How- 
ever, this is essentially a mean-field argument, which as- 
sumes that the substrate and enzyme molecules are ran- 
domly distributed in space at all times. Fig. [6] shows 
that spatio-temporal correlations between the enzyme 
and substrate molecules can have a dramatic effect on 
the existence of bistability. When the enzymes' reacti- 



Figs. [3]|n] show that enzyme-substrate rebindings are 
significant when the concentration of enzyme and sub- 
strate is on the order of 100 nM, which is a biologically 
relevant range Fig. 84 oi the Supporting Infor- 

mation shows that when the concentrations of all species 
are increased by more than a factor 10 from those used in 
Fig. [5l the system becomes bistable. While the distribu- 
tion of rebinding times does not depend on the concen- 
tration, the competition between phosphatase and kinase 
molecules in the bulk for binding to the substrate does, 
in such a way that the system is driven deeper into the 
bistable regime (see Fig. S5 in Supporting Information). 
Increasing the concentration can thus overcome the effect 
of enzyme-substrate rebindings. 



IV. DISCUSSION 

Multi-site phosphorylation is omnipresent in biological 
systems. Perhaps the best known and arguably the most 
studied example is the dual phosphorylation cycle of the 
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FIG. 6: Relative Kpp concentrations as a function of 
the lifetime of the inactive state of the enzyme, Trd- 
[K] total = [Kpp] + [Kp] + [K] was increased to 500 nM to 
bring the system to a regime where the mean-field model 
based on the chemical rate equations predicts bistability. At 
each value of t^^i , the particle-based model was simulated un- 
til it reaches steady state, starting from four different initial 
conditions, [Kpp] /[K] total ~ (blue), 0.3 (green), 0.7 (red), 
and 1 (cyan). While the mean- field model shows bistability 
over the whole range of Ti-oi (black dotted lines) , the particle- 
based model exhibits a bifurcation from mono- to bistability 
at t « 100/xs. At this bifurcation point, the system critically 
slows down, as a result of which it does not even equilibrate 
after 350 s. 



MAPK pathway, studied here, but other well-known ex- 
amples are the Kai system [3a|, the CDK inhibitor Sicl 
[13, the NFAT system H, and the CAMKII system [s^. 
Multi-site phosphorylation can lead to an ultrasensitive 
respon se to a threshold response Q , to bistabil- 

ity p^l3g]7or synchronise oscillations of phosphorylation 
levels of individual protein molecules [36|, provided the 
enzymes act via a distributive mechanism. We have stud- 
ied using a particle-based model a dual phosphorylation 
cycle in which the enzymes act according to a distribu- 
tive mechanism. Our results show that rapid enzyme- 
substrate rebindings can effectively turn a distributive 
mechanism into a processive mechanism, leading to loss 
of ultrasensitivity and bistability. Moreover, our results 
reveal that enzyme-substrate rebindings can significantly 
speed up the response, with slower diffusion leading to a 
faster response. While rebindings have been predicted to 
affect the noise in signal detection [H, H^j, our results 
predict that they can also drastically change the macro- 
scopic behaviour of the system. 

Our results reveal that enzyme-substrate rebindings 
occur on short length and time scales. Rebindings are 
important up to time scales of about 1 — lOnis (Fig. 14]), 
corresponding to the time for a protein to diffuse over a 
few molecular diameters. Beyond those length and time 
scales the dissociated enzyme and substrate molecules 
have essentially lost memory where they came from, and 
they would have to find each other again at random. An 
important question is whether we should not have taken 
orientational diffusion into account, precisely because re- 



bindings occur at comparable length and time scales. 
However, the first and second phosphorylation site are 
often close to each other on the substrate, e.g. separated 
by only a single amino-acid residue j4l| . suggesting that 
enzyme-substrate rebindings can indeed occur without 
significant orientational diffusion. Moreover, our model 
does not include molecular crowding, and it seems likely 
that subdiffusion caused by crowding can significantly 
extend the time scale over which rebindings occur [42|. 

The importance of enzyme-substrate rebindings de- 
pends on the lifetime of the inactive state of the en- 
zymes. For a typical protein diffusion constant oi D = 
1 — 10 /im^s~^ jisL H. [isl] ■ the rebinding probability drops 
below 10 % when the enzyme re-activation time becomes 
longer than 10ms (Fig. HP). Slow enzyme re-activation 
may thus be critical for generating bistability and ultra- 
sensitivity. To our knowledge, re-activation times of en- 
zymes in MAPK pathways have not been measured yet. 
The re-activation time of a kinase will depend sensitively 
on the order in which ADP and modified substrate dis- 
sociate from it, and ATP and substrate bind to it. If 
a kinase can bind its substrate irrespective of the nu- 
cleotide binding state, then nucleotide exchange will not 
be rate limiting. If the ADP is released before the modi- 
fied substrate, but ATP binding is required for binding of 
the next substrate, then ATP binding might be the rate- 
limiting step; with mM ATP concentrations, this is how- 
ever expected to yield fast re-activation times, of order 
microseconds. If the modified substrate must dissociate 
before ADP can dissociate and ADP must dissociate be- 
fore the kinase can bind substrate again, then the rate of 
ADP release may become rate limiting. A recent study 
on a protein kinase provides support for the latter sce- 
nario, with an ADP release rate that is on the order of 
100ms [321. This suggests that slow ADP release may 
allow for ultrasensitivity and bistability, although more 
work is needed to explore these mechanisms in depth. 
Concerning bistability, it is possible that bistability re- 
quires the phosphatase to act distributively [l^. Bista- 
bility could thus be lost if the mechanism by which the 
phosphatase acts changes from a distributive to a pro- 
cessive mechanism due to rebindings. To our knowledge, 
it is unknown what the minimum time is to re-activate 
a phosphatase. It is conceivable that this time is very 
short. Rapid re- activation of the phosphatase could thus 
lead to loss of bistability. 

Our results show that experiments to determine 
whether an enzyme acts distributively or processively 
should be interpreted with care. These experiments are 
often performed by investigating the time courses of the 
concentrations of the intermediate and final products [1] . 
If the amount of intermediate products exceeds that of 
the enzyme, then the mechanism must be distributive. 
However, our results reveal that the converse does not 
necessarily imply that the mechanism is processive, as 
commonly assumed: enzyme-substrate rebindings can 
turn a distributive mechanism into a processive one, with 
the concentration of the intermediate product remain- 
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ing below that of the enzyme. We stress that the ques- 
tion whether an enzyme acts processively because of re- 
bindings or because it remains physically attached to the 
substrate is biologically relevant, because the importance 
of enzyme-substrate rebindings strongly depends on the 
conditions. It depends on the diffusion constants of the 
components, the lifetime of the inactive state of the en- 
zyme, and on the concentrations of the components. All 
these factors may vary from one place in the cell to an- 
other and will vary from one cell to the next. In fact, an 
enzyme that operates according to a distributive mecha- 
nism in the test-tube may act processively in the crowded 
environment of the cell. 

Finally, how could our predictions be tested experi- 
mentally? If the enzyme of interest is a kinase, then one 
experiment would be to change the lifetime of the inac- 
tive state by varying the ATP concentration or by mak- 
ing mutations that change the ADP release rate. An- 
other proposal would be to study the enzyme kinetics 
as a function of the concentration of a crowding agent, 
such as PEG [S^. Crowding will slow down diffusion, 
and will, because of subdiffusion [i^l, increase the time 
that an enzyme and a substrate molecule that are in close 
proximity, stay together. Both effects will make enzyme- 
substrate rebindings more likely. Studying the input- 
output relation and the time course of the intermediate 
and final products [1, [s^] for different levels of macro- 
molecular crowding will shed light on the importance of 
spatio-temporal correlations for the macroscopic behav- 
ior of biological systems employing multi-site modifica- 
tions. 



V. METHODS 

A. Green's Function Reaction Dynamics 

A reaction-diffusion system is a many-body problem 
that can not be solved analytically. The key idea of 
GFRD is to decompose the many-body problem into sin- 
gle and two-body problems, which can be solved ana- 
lytically using Green's functions [H, . These Green's 
functions are then used to set up an event-driven algo- 
rithm, which makes it possible to make large jumps in 
time and space when the particles are far apart from 
each other. In the original version of the algorithm, the 
many-body problem was solved by determining at each 
iteration of the simulation a maximum time step such 
that each particle could interact with at most one other 
particle during that time step [l^, ■ In the enhanced 
version of the algorithm presented here, called eGFRD, 
spherical protective domains are put around single and 



pairs of particles [4g|. This allows for an exact, asyn- 
chronous event-driven algorithm (see Supporting Infor- 
mation) . 

B. MAPK model 

The model of the distributive, MAP kinase dual phos- 
phorylation cycle is sketched in Fig. [l]and described by 
Eqs. [Uni The rate constants are fci = 0.027 nM^^ • s"\ 
k2 = 1.35s"S h = 1.5s-\ ki = 0.056 nM-i • s-i, 
ks = 1.73 s'l, fee = 15.0 s-i, kj = In2/Trei. The 
protein diameter a = 5 nm. ki and are the in- 
trinsic association rates, which are the association rates 
for two species in contact; fc2 and k^ are the intrin- 
sic dissociation rates While in the particle-based 
model the diffusion of the particles is simulated explic- 
itly, in the mean-field model based on the ODE chem- 
ical rate equations, diffusion is described implicitly by 
renormalizing the association and dissociation rates [33| : 
l/kon = 1/fca-l-l/fcD and l/kos = l/^d + -fi^eq/te, where 
kon and fcoff are the renormalized association and dissoci- 
ation rates, respectively, k^^ = ki, and k^ = k2, k^ are 
the respective intrinsic association and dissociation rates, 
/cd = 4:naD is the diffusion-limited association rate, and 
Kcq = k^/kd = kon/koS is the equilibrium constant. The 
particles were put in a cubic volume of 1 /im"^ with pe- 
riodic boundary conditions. The total enzyme concen- 
tration [KK] -I- [P] is 100 nM corresponding to 60 copies 
of molecules in the volume, and the total substrate con- 
centration [K] -I- [Kp] + [Kpp] is 200 nM or 120 copies of 
molecules in Figs 3, 4 and 5, and 500 nM or 300 copies 
of molecules in Fig 6. The processive model consists of 
the following six reactions, sharing the same rate con- 
stants as the distributive model: KK -I- K ^-^^ KK — K ^ 
KK-Kp KK+Kpp, P+Kpp ^ P-Kpp P-Kp ^ 
P + K. 
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1 Enhanced Green's Function Reaction Dynamics 



1.1 Overview 

We present the enhanced Green's Function Reaction Dynamics (eGFRD) simulation algorithm. 
We provide the concepts required to understand the outline of the algorithm, but details on the 
algorithm, such as the actual mathematical expression for the employed Green's functions, other 
numerical procedures, and performance analyses, will be given in a forthcoming publication [T]. 

We solve the many-body reaction-diffusion problem by decomposing it into a set of many one 
body (single) and two body (pair) problems, for which analytical solutions (Green's functions) exist. 
In the original version of the GFRD algorithm, the many-body problem was solved by determining 
at each step of the simulation a maximum time step At such that each particle could interact 
with at most one other particle during that time step. In practice, the maximum time step At 
was determined as follows: 1) an interaction sphere of radius H^/GDjAi is drawn around each 
particle, where Di is the diffusion constant of the i-th particle, and H is a user-set error control 
parameter (usually 3 or 4); 2) the maximum time step At is then set by the requirement that 
each interaction sphere overlaps at most with one other interaction sphere. Subsequently, for each 
single particle and pair of particles a tentative reaction time is drawn, after which all particles are 
propagated simultaneously up to the smallest tentative reaction time, or to the maximum time 
step if that is smaller than the smallest tentative reaction time [21 [3] . Although already up to five 
orders more efficient than conventional reaction Brownian Dynamics [3] and also very accurate by 
its own right, the original GFRD algorithm has three major drawbacks; 1) due to the synchronous 
nature, the decomposition into one and two-body problems has to happen at every simulation 
step; 2) all components in the system are propagated according to the smallest tentative reaction 
time, making the performance sub-optimum; 3) the decomposition into single particles and pairs 
of particles involves cut-off distances, which makes the algorithm inexact. The systematic error 
is controlled by the H parameter, which determines the probability that during a time step At a 
particle travels a distance further than the maximum distance set by the requirement that each 
particle can interact with at most one other particle. This means that there is a trade-off between 
performance and error. 

In the current work, we overcome the drawbacks of the original GFRD scheme by putting 
protective domains around single particles and pairs of particles [5]. In this scheme, the next event 
of a domain can either be a reaction, or one particle leaving the domain. The tentative exit time 
for the latter event is computed by imposing an absorbing boundary condition on the surface of the 
protective domain. This makes the algorithm exact, and allows for an asynchronous event-driven 
algorithm. 

In the following sections, we explain how the reaction-diffusion problem in a spherical protective 
domain is solved for the one-body (Single problem) and two-body case (Pair problem), and then 
describe how the simulations of the different domains are integrated. 

1.2 Single particle events 

We consider a single particle of diameter d surrounded by a spherical protective shell of radius a 
(Fig. ISll fA)). Motion of a freely diffusing spherical particle is described by the Einstein diffusion 
equation, 

dtpi{r,t\ro,to) = DVV(r, t|ro, to), (1) 
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Figure SI: Single and Pair objects. To solve the many-body problem exactly, protective domains 
are put around single particles (A) and pairs of particles (B). (A) The radius of the protective 
domain of a single particle is denoted by a. (B) To solve the reaction-diffusion problem of two 
particles that can react with each other and diffuse within a protective domain with radius TZ, we 
construct two protective domains: one for the center-of-mass R, with radius or, and one for the 
inter-particle vector r, with radius Ur- The radii an and ar can be freely chosen, provided that 
when R and r would reach their maximum lengths, i.e. when |R| = an and |r| = a^, the particles A 
and B would remain within the protective domain. The latter means that an and ar should satisfy 
the following two constraints: 1) afi + arD^/ {Da + Db) < ^ — which reflects that particle A 
should remain with the protective domain with radius TZ; 2) aR + arDs / [Da + Db) < 1Z — (Tb/2, 
reflecting that B should remain with the protective domain with radius TZ. Although qr and 
ar can be freely chosen provided that these constraints are met, an efficient choice is given by 
a^/-D|j = (ar — ro)^/-D,^, meaning that the average time for R to reach the boundary of its domain 
by free diffusion, equals that of r. To illustrate the constraints, panel B shows a scenario where R 
and r reach their maximum lengths; here. Da = Db- 



where the Green's function pi(r, t|ro, to) denotes the probability that the particle is at position r 
at time t given that it was at rg at time to- We obtain the Green's function pi{r, t|ro, to) by solving 
the diffusion equation, Eq. [H with the following initial and boundary conditions 

pi(r,to|ro,to) = 5(r-ro), (2) 
pi(|r - rol = a,t|ro,to) = 0, (3) 

where 5 denotes the Dirac delta function. 

From the Green's function one can obtain the survival probability 

Si{t\ro,to) = / drpi(r,t|ro,to), (4) 

J |r-ro|<a 

which is the probability at time t that the particle remains within the protective sphere of radius 
a. This is related to the probability per unit time that the particle escapes the domain for the first 
time, 

(?™(t|ro,to) = -55(t|ro,to)/9t. (5) 

Sampling from this escape-propensity function yields a tentative escape time tescape- 

It is also possible that the particle undergoes a unimolecular reaction. The probability that the 
next reaction happens in an infinitesimal time interval t and t + dt is [6] 

^rcaction^^l^^^^^ = A;exp(-A:(t - to))dt, (6) 
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where k is the first-order reaction rate. This distribution can be used to obtain the next tentative 
reaction time treaction- 

The next event time of a Single is given by the smallest of the two tentative event times, namely, 



^single — mill(^escape) ^reaction)' (7) 



1.3 Particle pair events 



To describe the diffusion and the reaction of a pair of particles, we use the distribution function 
P2(ryi) rs, i|r^O! r_BO) *o)) which gives the probability that the particles A and B are at positions 
and at time t, given that they were at r^o and r^o at time to. This distribution function 
satisfies for |r| > a, where a = {dA + d_B)/2 is the cross-section with cIa and cIb the diameters of 
particles A and B, respectively, the following diffusion equation: 

dtP2{rA,rB,t\rAo,rBo,to) = [DaV\ + DbV^] P2(rA, r^, tlr^o, rso, ^o)- (8) 

We aim to solve this equation for two particles that can react with each other and diffuse within a 
protective domain. To our knowledge, it is impossible to solve this equation and obtain the Green's 
function directly. We therefore apply the following tricks. 
First, we make a coordinate transformation 

r = tb - rA, (10) 

and define the operators 

Vr = d/dR, (11) 

Vr = d/dr. (12) 

Eq. [8] can then be rewritten as: 

9tP2(R,r,t|Ro,ro,to) = [i?i?V|i + I^rV?] p2(R, r, t|Ro, tq, to), (13) 

where = DaDb/{Da + Db) and Dr = Da + Db- This equation describes two indepen- 
dent random processes, one for the inter-particle vector r and another for the center-of-mass vec- 
tor R. This means that the distribution function ^2(1"^; r^, ^Ir/io, rso, *o) can be factorized as 
p^{'R,t\Ilo,to)p2{T^,t\ro,to) and that the above equation can be reduced to one diffusion equation 
for the coordinate R and another for the coordinate r: 

dtpf{K,t\Ko,to) = DRVlj>f(R,t\Ko,to), (14) 
dtpl{r,t\ro,to) = DrVlpl{r,t\r,to). (15) 

The crux is now to define one protective domain for the interparticle vector r, with radius a^, 
and another for the center-of-mass vector R, with radius ur (see Fig. ISl( B)). These domains have 
to be chosen such that when the inter-particle vector r and the center-of-mass vector R would 
reach their maximum lengths, given by |r| = a,, and |R| = api, respectively, the particles A and B 
would still be within the protective domain for the two particles. 
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The diffusion equation for the center-of-mass vector now has to be solved with the boundary 
conditions 

pfiK,to\Ko,to) = <5(R-Ro), (16) 
p^(|R-Ro| =ai?,t|Ro,to) = 0. (17) 

This problem, of the center-of-mass diffusing in its protective domain, is similar to that of the 
single particle diffusing in a protective domain as discussed in the previous section. From the 
corresponding propensity function g2^(t|ro) we can draw a tentative time t^j at which the center- 
of-mass leaves its protective domain. 

The solution for the diffusion equation for the inter-particle vector is less trivial, since it should 
take into account not only that the inter-particle vector can leave its protective domain, but also 
that the two particles can react with each other. This reaction is included as an extra boundary 
condition, yielding the following boundary conditions for the inter-particle vector r: 

5(r-ro), (18) 
0, (19) 

kaPli\r\=a,t\ro,to), (20) 

where d/dr denotes a derivative with respect to the inter-particle separation r. Eq. [20] is the 
boundary condition that describes the possibility that A and B can react with a rate ka once they 
are in contact. Here, j(cj, t|ro, to) is the radial flux of probability ^2(1^1 t\^o,to) through the "contact" 
surface of area Aira'^. This boundary condition, also known as a radiation boundary condition [7], 
states that this radial flux of probability equals the intrinsic rate constant ka times the probability 
that the particles A and B are in contact. In the limit /cq — > 00, the radiation boundary condition 
reduces to an absorbing boundary condition P2(kl ~ <^i^|roiio) = Oi while in the limit ka the 
radiation boundary condition reduces to a reflecting boundary condition. 

From the Green's function for the inter-particle vector r, p2{r,t\ro,to), we can obtain two 
important quantities. The first is the time tbimo at which the inter-particle vector crosses the 
reaction surface given by |r| = a — meaning that the particles A and B react with each other — and 
the other is the time tr at which it "escapes" through the boundary of the protective domain given 
by |r| = Ur- The time at which the next event happens, be it a reaction or an escape, can be 
obtained through the survival probability, which is given by 

55(t|ro,to) = / drpl{r,t\ro,to). (21) 

J <T<|r| <ar 

The propensity function g2(t|ro,to)) which is the probability that the next event happens between 
time t and t + dt, is related to the survival probability by 

r/,| , N _ dS^{t\ro,to) 
q2it\ro,to) = — 

To know which of the two event types, reaction or escape, happens at time t, we split this quantity 
into two components, 

g^(t|ro,to)= qUt\ro,to)+q^^{t\ro,to) (23) 

= ijr|=a dSDrf^pl{r, t|ro, to) - ijr|=a, dSDrf^pl{r, t|ro, to), (24) 
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P2(r,io|ro,io) = 

pI{\t\ = a.r,t|ro,to) = 
d 

-j(cr,t|ro,to) = 47rcr^Dr o-P2(r,i|ro,to)||r|=f7 = 



(22) 



where in the first term dS denotes an integral over the reaction surface at |r| = a, and in the 
second term an integral over the boundary of the protective domain |r| = Ur- The reaction rate 
g2(i|rO)^o) is the probability that the next reaction for a pair of particles, initially separated by 
ro, occurs at time t, while the escape rate g2'^(f|ro, to) yields the probability that the inter-particle 
distance reaches and escapes from the protective domain for the first time at time t. We can 
draw the tentative time t for the next event, be it an escape or a reaction event, from Eq. [22| and 
then determine which of the two takes place from the ratio of q2{t\rQ,tQ) and g2'^(t|ro, to) at time 
t. Alternatively, we can draw a tentative time for a bimolecular reaction, tbimoi from 52(^^0)^0) 
and a tentative time for an escape event, tr, from 52^^ (t|ro, to); which of the two events can occur 
is then the one with the smallest tentative time (see below). The function ■§pP2{^,'t\rQ,tQ) can be 
used to sample the exit points on the relevant surfaces. 

It is possible that the particles A and B do not only react with each other, but also can undergo 
a unimolecular reaction of the type X — > . . . . In the same way as in the Single problem (Eq. [6|) , we 
can also draw the times tmono,A and tmono,_B at which the particles A and B undergo a first-order 
reaction, respectively. 

The next event of a pair of particles in a protective domain is thus one of the following events: 1) 
the center-of-mass leaving its domain; 2) the inter-particle event leaving its domain; 3) a bimolecular 
reaction; 4) unimolecular reaction of molecule A; 5) a unimolecular reaction of molecule B. The 
event that actually takes place is the one with the smallest tentative time. The next event time for 
a protective domain with two particles is thus given by 

tpair — mill(tij , tf , tbiujoi tjnono,j4; tmonOjij) ■ (^^) 

1.4 Algorithm outline 

The outline of the eGFRD algorithm is given by: 

1. Initialize: Reset the simulator time (tsim ^ 0). For each particle in the system, draw a 
spherical protective domain of appropriate size. When two particles are very close, create a 
Pair between them. Otherwise, create a Single object for each of the particles. Then, for 
each of the Single and the Pair objects, draw the next event type and the next event time 
according to the formulations in the previous sections, and chronologically order the events 
in the scheduler. 

2. Step: Pick the next event with the smallest scheduled time t from the scheduler. Update the 
simulator time tsim t. 

• Single event 

— If the event is a Single escape event, then (1) propagate the particle to a randomly 
determined exit point on the surface of the protective domain; (2) check if there are 
protective domains that are close to the new position of the particle; (3) if there are, 
burst the neighboring domains, and propagate the particles in the burst domains 
to a new position, and check if the current Single particle can form a Pair with 
one of the neighboring particles; (4) if a Pair is formed, discard the current Single, 
determine the new Pair event time (Eq. [25|) . and schedule the new Pair event on 
the scheduler; (5) for each of the particles contained in the stepping Single or the 
burst domains that are not used in formation of the Pair, draw a new domain and 
schedule a Single event on the scheduler. 
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— If the event type is Single reaction, (1) propagate the particle to a point r within 
the protective domain according to the Single Green's function pi(r, isim|riast) ^last), 
where tiast is the time the Single was created or the last time it stepped, and riast is 
the position of the particle at tiasti (2) execute the reaction by replacing the particle 
with one or more of the product particles placed next to each other; (3) for each 
of the newly created particles, draw a new protective domain and schedule a Single 
event on the scheduler. 

• Pair event 

— If the event type is Pair reaction, meaning that the two particles in the domain 
react, (1) draw the new R from p^(R, tsimlR-iasti ^last), where Riast is the position of 
R at tiast) which is the time at which the Pair was formed; (2) remove particles A 
and B of the Pair from the simulator; (3) place the product particle(s) at the new 
position R; (4) draw protective domain(s) around the new particle(s), and schedule 
Single event (s) on the scheduler. 

— If the event type was r escape, meaning that the inter-article vector r leaves its 
protective domain, (1) sample the new R position as above; (2) sample the r exit 
point from ^^2' (^) determine the new positions of A and B, ta and re, by putting 
the R as calculated in (1), and the exit point r on the surface of the inter-particle 
protective domain as calculated in (2), into Eqs [9l and \W\ (4) delete the Pair, (5) 
create a Single domain and schedule a new Single event on the scheduler for both A 
and B. 

— If the event type is R escape, meaning that the center-of-mass leaves its domain, 

(1) sample the new inter-particle vector r with P2(''' ^sim|riast) *iast)5 where rjast is 
the inter-particle vector at the time tiast it was last updated; (2) sample the R exit 
point from ^p^; (3) displace the particles A and B to the new positions; (4) delete 
the Pair; (5) create two Singles and schedule them on the scheduler. 

— If the event type was single reaction, (1) burst the pair domain and update the posi- 
tions of the particles A and B by sampling p^(R, tsim|R-iast, ^last) andp^(r, tsim|riast, ^last); 

(2) execute the reaction of the reacting particle according to the same procedures as 
used in the Single event; (3) create a new Single domain for the other (non-reacting) 
particle and schedule it on the scheduler. 

3. Go to 2. 

It can happen that more than two particles come very close to each other, making it difficult 
to draw protective domains of sufficient size around each of the particles [5]; this could bring the 
simulations to a standstill. To preempt this scenario, the algorithm puts one protective domain 
around the "squeezed" particles to form a third type of object called Multi. The particles in 
this domain are propagated according to Brownian dynamics [4] until the particles recover from 
the squeezed condition. Since it is guaranteed that Brownian dynamics converges to the correct 
solution when a sufficiently small step size is used[4j, this squeezing recovery procedure does not 
affect the overall accuracy of the simulation. 

The actual forms of the single and pair Green's functions, efficient numerical evaluation methods 
for the Green's functions, more details on the algorithm including the recovery procedure from 
squeezing, and handling of surfaces will be described in a forthcoming publication [1]. 
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2 Rebinding-time distribution 



In this section we present scaling relations for the rebinding-time distributions of two particles that 
can diffuse and react with each other in a large compartment. These results give a mathematical 
interpretation of the non-monotonic form of the enzyme-substrate rebinding-time distributions, 
shown in Fig. 4 (panels A and B) of the main text. 

The problem is reduced to solving the reaction-diffusion equation for a random walker that can 
diffuse in a domain internally bounded by a sphere of radius u, and to which it can bind with an 
intrinsic rate /ca once it is in contact with the sphere. This represents the evolution of the inter- 
particle vector r describing the distance between a substrate molecule and the enzyme molecule 
that has just modified it. 

The rebinding probability can be obtained from the Green's function P2(rj i|ro, io)- The proba- 
bility that a particle that starts at the origin given by |r| = a returns to the origin at a later time 
t, is given by 

(J _ gi3t/a2{i+M'^^(i + ha)ex{c ((1 + ha) 



' Dt 

= wv/io; ^ —' P"' 

where erfc is the complementary error function. If we assume that the enzyme becomes active 
immediately after enzyme-substrate dissociation, then, according to Eq. [20} the rebinding-time 
probability distribution can be expressed as 

p'^.^it) = k,p2{a,t\a,0). (27) 

This rebinding-time distribution has a number of properties. Firstly, the total probability that 
there is a rebinding is smaller than one: 

° 1 

dtPrehit) = 4,^, < 1. (28) 

2 

Secondly, upon a variable change t = r— — r"77T— =--ttt , the rebinding probability distribution can 
be rescaled as 

P^ehir) = Y^^j^fir), fir) 



e^erfc (\/t) 



(29) 

The function /(r) has the shape of two power laws /(r) ~ for r <C 1 and /(r) ~ ^^^^ for r ^ 1, 
in accordance with the results presented in Figure 2 of the main text. In fact, we can estimate 

2 

the time for the inflection point to be Tmoi = tt, — , . ^ Here, Tmoi represents the time after 
which most of the rebindings correspond to particles that start at contact, but wander away from 
the reaction sphere before they return to and rebind the reaction sphere. We stress that these 
trajectories are rebinding trajectories: we thus exclude trajectories where particles diffuse in the 
bulk and come back in a memory-less fashion (see also below). The scaling for t < t^oI can 

be understood by noticing that on this time scale particles stay close to the surface of the reaction 
sphere; indeed, on this time scale the particles essentially see a flat reaction surface, meaning that 
the return-time distribution is that of a ID random walker as described in the main text. The 
scaling for t > Tmoi can be understood by observing that on this time scale the particles 
have diffused away from the surface of the sphere; the particles now see the entire sphere, which 
means that the rebinding-time distribution is that of a 3D random walker returning to the origin. 
Interestingly, Tmoi depends on fca- When is increased, the probability that a particle binds 
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the target upon contact increases. Hence, the probabihty that after a time t the particle is still 
performing a ID random walk close to the surface, decreases as fca increases — the particle has either 
reacted with the surface, or escaped from the surface, thus performing a 3D random walk. 

So far we have assumed that upon dissociation, the enzyme and substrate can rebind as soon 
as they are in contact again. If, however, they can only rebind after the enzyme has become active 
again, then the rebinding-time distribution is given by 

Preh{t)= / dt' / i7Tr^drp*{r,t'\a,0)k^p^^tit')p2iT,t-t'\r,0), (30) 

Jo Ja 

where 

Pact(i) = A^acte-^-** (31) 

is the enzyme reactivation distribution with /cact = l/''"reh and P2{r,t\rQ,Q) is the solution of the 
Smoluchowski equation with a reflecting boundary condition — this reflects the idea that when the 
enzyme has not become active yet, the substrate cannot bind it. 

While we do not know an analytical expression for the rebinding-time distribution of Eq. [30l we 
can derive a lower bound p™™(t) and an upper bound p™^(i) for it, such that p™^{^{t) < Preh{t) < 
p™g'^(t). The upper bound p^^^{t) is based on the inequality 

P2{a, t - t'\r, 0) < p2{(7, t - t'\a, 0), (32) 

with equality for r = a. This inequality expresses the fact that the probability that the particles 
are in contact at a later time t decreases with the initial distance. This yields the upper bound 

PT^{t) = fdt' k,p^t{t')p2{<T, t - t'\a, 0). (33) 
Jo 

Using the solution for p2{(J,t — t'|cr, 0) as described in [29] one can show that for small t the bound 
increases with t as while for large t it decreases with t as The upper bound p™b^(t) 

thus has a non-monotonic behavior, going to zero at both zero and infinity. The position of the 
maximum depends on the two time scales Trei = l/^act and t^o\ and is located in the interval 
[MIN(r„oi, r,ei) , MAX(r„,oi, r^ei)] • 

The lower bound p^^{^{t) is based on the inequality P2{f'-, t\a, 0) < P2{r, t\a, 0). This reflects the 
idea that with a radiation boundary condition, the particle can react with the reactive sphere (and 
thus leak out of the system), while with a reflecting boundary condition it cannot. This yields the 
following expression for the lower bound 

Cb(i) = fdt'Kp,,,{t')p2{a,t\a,0), (34) 
Jo 

= (l-e-'=-**)te(a,tk,0), (35) 

= 5act(t)p?eb(t), (36) 

where Sact(i) is the probability that the enzyme is active after time t and p^ebC^) enzyme- 
substrate rebinding-time distribution assuming that the enzyme is active at all times. This lower 
bound has a number of interesting scaling regimes, depending on the relative values of r^ei and Tmoi- 
They can be understood intuitively by making the following observations: as discussed above, for 
t < Tjnoi, Preb(^) ~ t~'^/'^, while for t > Tmoi, P^chi^) ~ t'^^'^] moreover, for t < Trd, Sact(t) ~ t, 
while for t > Trei, Sact(i) 1- Hence, for t < MIN(rmoi, Trei), P^^{t) ~ t'^/^ ^ ^ ^+1/2^ ^j^j^^ f^^. 
t > MAX(rn,oi, Trei), ~ t-^^^- Moreover, when r^ei < r^oi, P^^{t) ~ t'^^^ for r^ei < t < t^^u 
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Figure S2: The upper bound Eq. [33] (dashed hues) and lower bound Eq. [36] (sohd hues) of the 
rebinding-time distribution, given by Eq. [30] for three different scenarios: 1) Trd ^ t^oI (blue 
lines); Trei ~ t^oI (green lines); T^ei ^ "^moi (red lines). It is seen that both bounds converge when 
t > Trel- It is also seen that the difference between the bounds is rather small when t < Tmoi' 

The 

difference between the upper and lower bounds arises for t^^i < t < Tj-ei, when t^^i S> Tmoi (red 
lines). This is because in this case the reactive sphere (i.e., enzyme) is mostly still inactive, and 
the (substrate) particle thus tends to diffuse away from it. This phenomenon is captured by the 
lower bound, but not by the upper bound. For a comparison with the simulation results, see Fig. 

m 

because Sact — 1 and Preb(^) ~ t"^/"^. And in the scenario that Tmoi < Tre\, P'^hi^) ^ t^^^"^ for 

Tmol < t < Trel, bccaUSC Sact ~ t and P°eb(*) ~ t'^^'^ ■ 

Fig. IS2I shows the predictions of the upper bound Eq. [33] and lower bound Eq. [36] for the 
rebinding-time distribution given by Eq. [30], for three different scenarios: 1) Tj-d <^ t^qi (blue line); 
2) Trel ~ Tjnoi (green line); 3) t^cI > ''"mol (red line). It is seen that in all 3 scenarios the upper 
and lower bound for Preh{t) converge for t > T^ei- Indeed, in this regime, where the enzyme is 
active, Preh{t) scales as t~^/'^. It is also observed that when Tj-ei < Tmoi (blue and green lines), the 
difference between the upper and lower bound for Prch{t) is very small, even when t < Trei. This 
implies that both bounds are good approximations for Prcb(Oj we can thus conclude that, to a 
good approximation, Preb(i) scales as t^^"^ for t < MIN(trei, Tmoi) and t~^^'^ for Trei < t < t^oI when 
Trel < Tmol- The difference between the bounds arises when Tj-ei > t^^i (red line); in this scenario, 
the bounds differ in the regime Tmol < t < Trei- The question arises which bound is closer to the 
actual rebinding-time distribution, Preh{t)- To this end, we compare our analytical results with the 
simulation data, shown in Fig. IS3I 

In Fig. [S3] we show the simulation results of Fig. 4 of the main text. When t > Tbuik the 
rebinding-time distribution is exponential. This is due to particles that come from the bulk, and 
bind the reactive sphere in a memory-less fashion. This regime is not described by the analysis 
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Figure S3: The enzyme-substrate association-time distribution of Fig. 4 of the main text, to- 
gether with the scahng regimes as predicted by the analysis of the upper and lower bounds for the 
rebinding-time distribution (see Fig. [S2]) : in panel A Tj-ei « t^oI, while in panel B Tj-ei ^ t^oI- For 
t > Thuik, the association-time distribution is exponential, because on this time scale the particles 
meet each other at random in the bulk. As predicted by the analysis of the upper and lower 
bounds for the rebinding-time distribution (see Fig. IS2p . the enzyme-substrate association-time 
distribution scales as t^^'^ for t < MIN(Tinoi, t^c\), and as for MAX(rmoij "^rei) < t < Tbuik! while 
the t^/^ scaling is seen in both panels, the is only seen in panel A, because in panel B T\^u\k 

approaches Tj-ei. Panel B shows that the lower bound Eq. [36] correctly predicts the t~^^'^ scaling for 

Tmol < t < Treh when T^l > T^oh 

discussed above, which is performed for the geometry of an infinite, spherical domain internally 
bounded by a reactive sphere. For this geometry, in three dimensions, there is a probability that 
the particle escapes to infinity without returning to and reacting with the reactive sphere. In a 
bounded domain, when a particle escapes from the vicinity of the reactive sphere into the bulk, it 
will return to the reactive sphere on a time scale Tbuik- Therefore, the rebinding-time distributions 
for the infinite domain analyzed above and the finite domain of the simulations, are the same, 
but only up to Tbuik- This also means that to observe the different power-law scaling behaviours, 
Tbulk > MAX(rrei ) ^mol ) • 

In panel A of Fig. [S3l Tj-d ~ Tmol) while in panel B r^ei ^ Tmol- In both scenarios Preh{t) ~ (t)^^'^ 
when t < MIN(rj.ei, Tjnoi), in accordance with the analysis of the upper and lower bounds of Preb(^) 
presented above. Both panels also show that when MAX(rjnoi, Trd) < t < Tbuik) Preb(*) ~ t~^/'^. 
An interesting regime is Tmol < t < Trei in the case that Trei > Tmoi (panel B). It is seen that the 
simulation results suggest that Preh{t) ~ t"^^"^ in this regime. This is predicted by the lower bound 
for Preh{t), Eq. [361 but not by the upper bound, Eq. [33] (see Fig. [S2l) . This can be understood by 
noting that in this regime, Tmoi < t < Tj-ei, the enzyme is mostly still inactive and the particle can 
thus diffuse away from the reactive sphere; while the lower bound of Eq. [36] captures this effect, the 
upper bound of Eq. [33|does not, since it is based on the inequality p2{cr,t — t'\r, 0) < p2{cr,t — t'\a, 0). 
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3 The effect of concentration 



Fig. IS4I shows the input-output relation as a function of concentration. Here, the concentrations 
of all components are increased by the same factor from the base-Hne values used in Fig. 5 A of the 
main text. It is seen that both the particle-based model and the mean-field model predict that an 
increase in concentration induces bistability, although the concentration at which the bifurcation 
occurs is higher in the particle-based model. 




Figure S4: Steady-state input-output relations for different concentrations. The concentrations of 
all components are increased by the same factor. (A) Baseline parameter values; the concentrations 
equal those corresponding to Fig. 5(A) in the main text: [K]totai = 200nM, [KK] -|- [P] = lOOnM). 
(B) 3x concentration ([K]totai = 600nM, [KK] + [P] = 300nM), (C) lOx concentration ([Kjtotai = 
2iiM, [KK] + [P] = l^M), (D) lOOx concentration ([K]totai = 20nM, [KK] + [P] = WfiM). 

Fig. IS5I elucidates the origin of why an increase in concentration can induce bistability, both in 
the mean-field model and the particle-based model. Bistability arises when a substrate molecule 
that has been phosphorylated once, is more likely to be dephosphorylated again than to become 
fully phosphorylated (similarly, the probability that after a fully phosphorylated molecule has 
been dephosphorylated once becomes fully phosphorylated again, should be higher than that it 
becomes fully dephosphorylated). We therefore plot in Fig. [S5] the probability that a substrate 
that has just been phosphorylated once, either binds the same kinase molecule as the one that just 
phosphorylated it (this is most likely due to a rebinding event), another kinase molecule (from the 
bulk), or a phosphatase molecule (from the bulk); the system is in a state where most substrate 
molecules are unphosphorylated. It is seen that the fraction of rebindings is fairly constant. This 
can be understood as follows: 1) the probability that a molecule returns to the origin before it 
looses memory where it came from is independent of the concentration (see Fig. 2 of the main 
text) — only the memory-less returns from the bulk depend on concentration; 2) when a rebinding 
event happens, it happens very fast: as Fig. 2 of the main text shows, rebindings are dominated 
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by events that occur on time scales of t < 1 ms. These time scales are so short, that the probability 
that an enzyme molecule from the bulk interferes with a rebinding event, is negligible, even up to 
concentrations of 100-1000 times the baseline value, i.e. 10 — 100 /xM; only above that concentration 
can molecules from the bulk effectively compete with those undergoing a rebinding trajectory, and 
will the probability that a dissociated molecule rebinds drop significantly. Up to a concentration 
of 10 — lOO/uM, there is thus an essentially constant probability, independent of the concentration, 
that both sites of a substrate molecule are modified by the same enzyme molecule. Now bistability 
can arise when the antagonistic enzyme in the bulk wins the competition from the agonistic enzyme 
undergoing the rebinding event and the other agonistic enzymes in the bulk. Fig. [S5] shows that 
when the concentration is increased, the competition between the kinase (the agonist) in the bulk 
and the phosphatase (antagonist) in the bulk changes in favor of the phosphatase. This is because 
the system is in a state where the substrate molecules are mostly unphosphorylated, and in this state 
the kinase molecules become increasingly sequestered by the unphosphorylated substrate molecules 
as the concentration is increased. This increases the probability that a molecule that has just been 
phosphorylated once, will bind a phosphatase (antagonist), which will drive it back towards the 
unphosphorylated state. Increasing the overall concentration thus changes the competition between 
the kinases and the phosphatases in the bulk in such a way that the driving force towards a state in 
which the substrate molecules are either fully unphosphorylated or fully phosphorylated, increases. 
In essence, increasing the concentration drives the system deeper into the bistable regime. This 
makes it possible to overcome the effect of rebindings, which tends to drive the system out of the 
bistable regime, as shown in Fig. 6 of the main text. 




Figure S5: The probability that a substrate molecule that has been phosphorylated once, will bind 
the same kinase molecule (blue) , another kinase molecule (red) or a phosphatase molecule (green) , 
for difference concentrations; the baseline values correspond to Fig.SA of the main text and Fig. IS4I 
It is seen that the fraction of events where the substrate molecule binds the same kinase molecule 
again is fairly constant, while the fraction of events in which the substrate molecule binds another 
kinase molecule strongly drops in favor of those in which the substrate molecule binds a phosphatase 
molecule, when the concentration is increased by a factor 10 from the baseline value — as shown in 
Fig. IS41 the system now becomes bistable. 
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